home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / testm.i < prev    next >
Text File  |  1996-02-12  |  6KB  |  215 lines
  1. /*
  2.    TESTM.I
  3.    Certify functions in fft.i and matrix.i
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func testm
  11. {
  12.   elapsed= old_elapsed= [0., 0., 0.];
  13.   write, "testing fft routines...";
  14.   timer, old_elapsed;
  15.   fft_test, 95;
  16.   fft_test, 32;
  17.   fft_test, 12;
  18.   timer, elapsed;
  19.   timer_print, "FFT elapsed time", elapsed-old_elapsed;
  20.   write, "testing LU decomposition routine...";
  21.   timer, old_elapsed;
  22.   testLU, 75;
  23.   testLU, 16;
  24.   timer, elapsed;
  25.   timer_print, "LU elapsed time", elapsed-old_elapsed;
  26.   write, "testing QR decomposition routine...";
  27.   timer, old_elapsed;
  28.   testQR, 75, 75;
  29.   testQR, 16, 16;
  30.   testQR, 25, 21;
  31.   testQR, 21, 25;
  32.   timer, elapsed;
  33.   timer_print, "QR elapsed time", elapsed-old_elapsed;
  34.   write, "testing SVD routine...";
  35.   timer, old_elapsed;
  36.   testSVD, 75, 75;
  37.   testSVD, 16, 16;
  38.   testSVD, 25, 21;
  39.   testSVD, 21, 25;
  40.   testSV, 21;
  41.   testSV, 64;
  42.   timer, elapsed;
  43.   timer_print, "SVD elapsed time", elapsed-old_elapsed;
  44. }
  45.  
  46. func fft_test(n)
  47. {
  48.   index= 2*pi*(indgen(n)-1.0)/n;
  49.   z= sin(index*3);
  50.   zf= fft(z, 1);
  51.   z3= z2= array(0i, n);
  52.   z3(4)= -0.5i*n;
  53.   z3(-2)= 0.5i*n;
  54.   zb= fft(z, -1);
  55.   if (max(abs(zf-z3))>1.e-12*n || max(abs(zb-conj(z3)))>1.e-12*n)
  56.     write, "***WARNING*** failed 1D fft test";
  57.   if (n<=96 || force2d) {
  58.     z*= cos(index*2)(-,);
  59.     zf= fft(z, [0, 1]);
  60.     z2(3)= z2(-1)= 0.5*n;
  61.     zb= fft(z, [], [-1, 0]);
  62.     if (max(abs(zf-sin(index*3)*z2(-,)))>1.e-12*n ||
  63.     max(abs(zb-conj(z3)*cos(index*2)(-,)))>1.e-12*n)
  64.       write, "***WARNING*** failed first 2D fft test";
  65.     zf= fft(z, 1);
  66.     zb= fft(z, -1);
  67.     if (max(abs(zf-z3*z2(-,)))>1.e-12*n ||
  68.     max(abs(zb-conj(z3)*z2(-,)))>1.e-12*n)
  69.       write, "***WARNING*** failed second 2D fft test";
  70.   }
  71. }
  72.  
  73. func TDcheck(c, d, e, b, x, s)
  74. {
  75.   check= _(   d(1)*x(1)    +    e(1)*x(2),
  76.        c(1:-1)*x(1:-2) + d(2:-1)*x(2:-1) + e(2:0)*x(3:0),
  77.                             c(0)*x(-1)   +   d(0)*x(0)   );
  78.   if (max(abs(check-b))>1.e-9*max(abs(b))) {
  79.     write, "***WARNING*** "+s+" tridiagonal solution doesn't check";
  80.     write, "   max relative error is "+pr1((max(abs(check-b)))/max(abs(b)));
  81.   }
  82. }
  83.  
  84. func testTD(n)
  85. {
  86.   c= random(n-1);
  87.   d= random(n);
  88.   e= random(n-1);
  89.   b= random(n);
  90.   TDcheck,c,d,e,b,TDsolve(c,d,e,b), "1D";
  91.   b2= random(n);
  92.   x= TDsolve(c,d,e,[b,b2])
  93.   TDcheck,c,d,e,b, x(,1), "2D(1)";
  94.   TDcheck,c,d,e,b2, x(,2), "2D(2)";
  95.   x= TDsolve(c,d,e,transpose([b,b2]), which=2)
  96.   TDcheck,c,d,e,b, x(1,), "2D(1)/which";
  97.   TDcheck,c,d,e,b2, x(2,), "2D(2)/which";
  98. }
  99.  
  100. func LUcheck(a, b, x, s)
  101. {
  102.   check= a(,+)*x(+);
  103.   if (max(abs(check-b))>1.e-9*max(abs(b))) {
  104.     write, "***WARNING*** "+s+" LUsolve solution doesn't check";
  105.     write, "   max relative error is "+pr1((max(abs(check-b)))/max(abs(b)));
  106.   }
  107. }
  108.  
  109. func testLU(n)
  110. {
  111.   a= random(n,n);
  112.   b= random(n);
  113.   LUcheck,a,b,LUsolve(a,b), "1D";
  114.   b2= random(n);
  115.   x= LUsolve(a,[b,b2])
  116.   LUcheck,a,b, x(,1), "2D(1)";
  117.   LUcheck,a,b2, x(,2), "2D(2)";
  118.   x= LUsolve(transpose(a),[b,b2])
  119.   LUcheck,transpose(a),b, x(,1), "t2D(1)";
  120.   LUcheck,transpose(a),b2, x(,2), "t2D(2)";
  121.   x= LUsolve(a,transpose([b,b2]), which=2)
  122.   LUcheck,a,b, x(1,), "2D(1)/which";
  123.   LUcheck,a,b2, x(2,), "2D(2)/which";
  124.   x= LUsolve(transpose(a),transpose([b,b2]), which=2)
  125.   LUcheck,transpose(a),b, x(1,), "t2D(1)/which";
  126.   LUcheck,transpose(a),b2, x(2,), "t2D(2)/which";
  127.   ai= LUsolve(a);
  128.   err= max(abs(ai(,+)*a(+,)-unit(n)))/max(abs(ai));
  129.   if (err>1.e-9) {
  130.     write, "***WARNING*** LUsolve inverse is fishy";
  131.     write, "   max relative error is "+pr1(err);
  132.   }
  133. }
  134.  
  135. func QRcheck(a, b, x, s)
  136. {
  137.   check= a(,+)*x(+);
  138.   err= max(abs(check-b))/max(abs(b));
  139.   if (err>1.e-9 && dimsof(a)(2)<=dimsof(a)(3)) {
  140.     write, "***WARNING*** "+s+" QRsolve solution doesn't check";
  141.     write, "   max relative error is "+pr1(err);
  142.   }
  143. }
  144.  
  145. func testQR(m,n)
  146. {
  147.   a= random(m,n);
  148.   b= random(m);
  149.   QRcheck,a,b,QRsolve(a,b), "1D";
  150.   b2= random(m);
  151.   x= QRsolve(a,[b,b2])
  152.   QRcheck,a,b, x(,1), "2D(1)";
  153.   QRcheck,a,b2, x(,2), "2D(2)";
  154.   x= QRsolve(a,transpose([b,b2]), which=2)
  155.   QRcheck,a,b, x(1,), "2D(1)/which";
  156.   QRcheck,a,b2, x(2,), "2D(2)/which";
  157. }
  158.  
  159. func SVcheck(a, b, x, s)
  160. {
  161.   check= a(,+)*x(+);
  162.   err= max(abs(check-b))/max(abs(b));
  163.   if (err>1.e-9) {
  164.     write, "***WARNING*** "+s+" SVsolve solution doesn't check";
  165.     write, "   max relative error is "+pr1(err);
  166.   }
  167. }
  168.  
  169. func testSVD(m, n)
  170. {
  171.   a= random(m,n);
  172.   s= SVdec(a, u, v);
  173.   if (anyof(s(dif)>0.0))
  174.     error, "***WARNING*** SVdec returned increasing singular values";
  175.   achk= u(,+) * (s*v)(+,);
  176.   sabs= max(abs(s));
  177.   err= max(abs(a-achk))/sabs;
  178.   if (err>1.e-9) {
  179.     write, "***WARNING***  SVdec decomposition doesn't check";
  180.     write, "   max relative error is "+pr1(err);
  181.   }
  182.  
  183.   s= SVdec(a, u, v, full=1);
  184.   if (anyof(s(dif)>0.0))
  185.     error, "***WARNING*** SVdec returned increasing singular values";
  186.   uu= u(,1:min(m,n));
  187.   vv= v(1:min(m,n),);
  188.   achk= uu(,+) * (s*vv)(+,);
  189.   sabs= max(abs(s));
  190.   err= max(abs(a-achk))/sabs;
  191.   if (err>1.e-9) {
  192.     write, "***WARNING***  SVdec decomposition doesn't check";
  193.     write, "   max relative error is "+pr1(err);
  194.   }
  195.   err= max(abs(u(,+)*u(,+)-unit(m)))+max(abs(v(,+)*v(,+)-unit(n)));
  196.   if (err>1.e-9) {
  197.     write, "***WARNING***  SVdec decomposition not orthogonal";
  198.     write, "   max relative error is "+pr1(err);
  199.   }
  200. }
  201.  
  202. func testSV(n)
  203. {
  204.   a= random(n,n);
  205.   b= random(n);
  206.   SVcheck,a,b,SVsolve(a,b), "1D";
  207.   b2= random(n);
  208.   x= SVsolve(a,[b,b2])
  209.   SVcheck,a,b, x(,1), "2D(1)";
  210.   SVcheck,a,b2, x(,2), "2D(2)";
  211.   x= SVsolve(a,transpose([b,b2]), which=2)
  212.   SVcheck,a,b, x(1,), "2D(1)/which";
  213.   SVcheck,a,b2, x(2,), "2D(2)/which";
  214. }
  215.